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ABSTRACT 


Considering  multifunctional  fibrous  composites,  laminates  and  structures,  a  unified, 
transformation  field  analysis  (TFA)  approach  is  developed  to  model  their  overall  response  to 
thermomechanical  loads,  as  well  as  local  phenomena  due  to  coupled  effects  and  damage.  Local 
stresses  and  strains  generated  as  a  result  of  piezoelectric  and/or  pyroelectric  behavior  of  the 
fibers  are  treated  as  transformation  fields,  which  remain  in  the  system  even  under  mechanical 
unloading.  Their  effect  on  the  stress  and  strain  distribution  at  the  microscale  of  a  unidirectional 
composite,  the  macroscale  of  a  laminate,  and  the  structural  scale  is  detennined  in  terms  of 
transfonnation  influence  factors  that  are  a  function  of  the  elastic  properties  of  the  fiber  and 
matrix  phases,  and  the  geometry  idealizations  at  the  three  length  scales.  Treatment  of  damage 
follows  the  same  scheme  but  the  transfonnation  fields  are  instead  determined  such  that  the  local 
stresses  in  the  affected  phase  are  removed.  The  result  is  tri-lengthscale  model,  which  seamlessly 
connects  the  stress  and  strain  fields  in  the  fiber  and  matrix  phases,  membrane  forces  and  bending 
moments  in  laminates,  and  the  overall  thennomechanical,  structural  loads  to  determine  the 
overall  response.  The  Mori-Tanaka  and  the  periodic  hexagonal  array  models  were  utilized  to 
model  the  micromechanical  behavior  of  the  individual  composite  laminas,  and  the  laminated 
plate  theory  was  applied  for  fibrous  laminates.  The  structural  response  is  determined  with  the 
ABAQUS  finite  element  code,  in  which  the  transformation  field  analysis  was  implemented  as  a 
user-defined  routine.  Examples  are  given  for  composite  laminates  and  structures  to  illustrate  the 
capabilities  of  the  developed  approach  in  predicting  the  overall  behavior  for  multifunctional 
composite  structures. 
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1.  INTRODUCTION 


In  the  rapidly  developing  field  of  multifunctional  composites,  robust  modeling  of  laminated 
structures  with  piezoelectric  fibers  is  lacking.  There  is  a  real  need  to  combine  analysis  of  large 
structural  parts  with  realistic  material  models,  which  represent  the  underlying  deformation 
phenomena.  There  is  no  shortage  of  work  in  multifunctional  composites  in  general,  but  the  focus 
has  been  mainly  on  finding  effective  properties,  including  mechanical,  electrical,  and 
electromechanical  constants.  Examples  are  found  in  the  work  of  Chen  (1994),  Aboudi  (2001), 
Berger  et  al.  (2005),  Kumar  and  Chakraborty  (2009),  and  Challagulla  and  Georgiades  (2011).  On 
the  other  hand,  a  combination  of  micro  and  macro  mechanical  modeling  of  laminated  plates  with 
embedded  piezoelectric  fiber  is  presented  by  Cook  and  Vel  (2012).  They  implemented  PFC 
actuators  composed  of  PZT-7A  fibers  embedded  in  an  epoxy  resin  inducing  shear  defonnation 
when  actuated  in  the  axial  direction.  Application  of  a  beam  bending  composed  of  a  PFC  shear 
actuator  embedded  between  two  graphite/polymer  substrates  is  found  in  their  work.  A  3D 
axisymmetric  laminated  shell  with  piezoelectric  sensors  and  actuators  is  presented  by  Santos  et 
al,  (2006).  The  shell  is  subject  to  mechanical  load  and  electrical  potential  on  the  external  surface. 
They  plotted  relations  between  stresses  and  displacements,  radial  and  axial,  and  the  applied 
electrical  potential  and  pressure.  Other  researchers  focused  on  the  implementation  of 
multifunctional  composite  laminates  in  structural  health  monitoring.  Diamanti  and  Soutis  (2010) 
used  PZT  transducers  to  act  as  actuators  and  sensors  to  detect  damage  in  a  carbon  fiber/epoxy 
laminate  due  to  low  velocity  impact.  Almost  the  same  technique  is  used  by  Giurgiutiu  and  Soutis 
(2011).  They  enabled  the  piezoelectric  wafer  active  sensors  (PWAS)  to  detect  damage  in  a 
composite  laminate  using  Famb  waves.  Huang  et  al,  (2009)  investigated  numerically  the  effect 
of  integrating  a  SHM  sensor  in  a  unidirectional  glass  fiber  laminate  on  the  onset  of  cracks  in  the 
resin  matrix.  The  same  problem  is  studies  experimentally  by  Ghezzo  et  al.  (2009).  Similarly, 
Fissenden  et  al.  (2009)  embedded  a  linear  array  of  piezoelectric  fibers  at  the  mid-plane  of  a 
carbon  fiber  reinforced  polymer  laminate.  They  simulated  actuation  of  the  PZT  fibers  by 
applying  a  radial  displacement  around  their  circumference  and  studied  the  effect  on  the  elastic 
strain  energy  density  of  the  system. 

This  work  focuses  on  the  micromechanical  analysis  of  a  structure  made  of  a  laminated 
composite  with  piezoelectric  fibers.  The  analysis  therefore  includes  the  interaction  between 
different  length  scales  of  the  problem.  From  the  broadest  to  the  finest,  there  is  a  structure  scale,  a 
laminate  scale,  a  composite  scale  and  a  phase  scale.  The  structure  scale  is  simulated  using 
ABAQUS  Explicit  fonnulation  while  the  laminate  scale  is  modeled  using  a  transformation  filed 
analysis  scheme.  The  composite  aggregate  is  solved  on  the  micromechanical  scale  using  a 
representative  volume  element,  which  models  a  periodic  hexagonal  distribution  of  the  fibers. 

This  report  describes  the  above  work,  which  was  completed  under  AFOSR  sponsorship.  The 
report  is  organized  as  follows:  Section  2  describes  the  tri-lengthscale  structural  analysis  model, 
which  combines  micromechanics,  a  laminate  theory,  and  the  finite  element  method.  Applications 
are  given  in  Section  3.  The  report  closes  with  conclusions  in  Section  4. 
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2.  A  TRI-LENGTHSCALE  MODEL  FOR  MULTIFUNCTIONAL  LAMINATED 

COMPOSITE  STRUCTURES 

Development  of  a  tri-lengthscale  model  for  multifunctional  composites  structures  was 
developed  in  two  steps.  First,  analysis  of  fibrous  laminates  with  an  underlying  micromechanical 
model  for  piezoelectric  laminas  was  developed  for  overall  membrane  forces  and  bending 
moments  using  a  Transformation  Field  Analysis  (TFA).  This  establishes  a  dual  lengthscale 
model,  which  relates  the  overall  response  of  a  fibrous  laminate  to  both  the  in-situ  properties,  and 
the  local  deformations.  Taking  this  analysis  to  the  next  level  involves  the  integration  of  the  TFA 
for  laminates  with  a  finite  element  framework.  The  model  is  described  in  detail  in  two 
publications  by  the  principal  investigator  (Bahei-El-Din  and  Micheal,  2012,  2013).  A  summary 
of  the  tri-lengthscale  model  with  relevant  equations  is  given  in  this  section. 

2.1.  MICRO-SCALE  MECHANICS  OF  A  FIBROUS  COMPOSITE 

Considering  a  unidirectional  lamina,  and  assuming  that  it  consists  of  Q  phases,  each  has  a 
volume  fraction  vr ,  the  phase  stresses,  er  ,  and  strains,  sr ,  are  related  to  their  overall 
counterparts  a,  £  ,  relative  to  the  lamina  coordinate  system  xk  ,  by: 

e=Y,°’£r-  (!) 

r=\,Q  r=l,Q 

Both  the  stresses  and  strains  appearing  in  Eq.  (1)  are  (6x1)  matrices  which  list  the  independent 
components  of  their  tensorial  counterparts.  Hence,  cr  =  [on,o22,o33,o2i,on,ol2],  and 

£  —  lAn ’£22 5 2^12 ]  • 

Constitutive  equations  of  the  phases  can  be  written  as  (Bahei-El-Din  and  Dvorak,  2000) 
crr  =  Lrer  +  Ar,  er  =  Mo  +  //,.  (2) 

Here,  Lr  and  Mr  =  Lr~x  are  elastic  stiffness  and  compliance  of  phase  r .  The  first  term  in  Eq.  (2) 
represents  the  response  to  mechanical  loads  and  the  second  term  represents  the  transformation  or 
eigen  fields,  //,.  and  Ar  =  -Lrjur ,  caused  by  effects  other  than  mechanical  loading.  Examples  of 

the  latter  are  thermal  changes  and  electrical  fields  induced  to  piezoelectric  phases.  The  phases  are 
assumed  to  be  transversely  isotropic.  In  terms  of  Hill’s  elastic  moduli  (Hill,  1964)  the  elastic 
stiffness  matrix  Lr  is  given  by 

k  +  m  k-m  i  0  0  0 
k  —  m  k  +  m  l  0  0  0 

t  in  0  0  0 

r~  0  0  0  /^00 

0  0  0  0  p  0 

0  0  0  0  0  m 


4 


In  tenns  of  engineering  moduli,  Hill’s  moduli  are  given  by: 
r  1  4  4v?l 

k  =  - - 1 - —  ,  £  =  2kv,  ,  n  =  EL+4kvL  =  EL+ 12  / k,  m  =  GT,  p  =  GL,  (4) 

Gj  ET  E  L 

where  EL,  Er  are  longitudinal  and  transverse  Young’s  moduli,  vL  is  Poisson’s  ratio 
corresponding  to  axial  straining,  and  GL ,  GT  are  longitudinal  and  transverse  shear  moduli. 

In  analogy  with  (2),  the  following  constitutive  relations  can  be  written  for  the  lamina  relative 
to  the  material  principal  axes  xk  , 

a  =  Ls  +  A,  J-Ma  +  Ji,  (5) 

where  L  and  M  =  L  are  overall  elastic  stiffness  and  compliance.  The  lamina  is  transversely 
isotropic  and  as  such  the  overall  elastic  stiffness  can  be  written  as  given  in  Eq.  (3)  in  terms  of 
overall  moduli  which  can  be  either  measured  or  modeled  (Bahei-El-Din  and  Dvorak,  2000). 

The  lamina  transfonnation  stress  A  and  strain  p  can  be  expressed  in  terms  of  their  local 
counterparts  as  (Dvorak  and  Benveniste,  1992) 

^  =  Z  VrA'rK-  P  =  X  UrBrJU^  (6) 

r= l,e  r=\,Q 

where  Ar  and  Br  are  strain  and  stress  concentration  factors,  which  provide  the  phase  fields  in 
tenns  of  their  overall  counterparts.  These  can  be  determined  from  a  material  model. 

Considering  overall  mechanical  loads  together  with  local  transfonnation  fields,  the  strains 
and  stresses  in  the  phases  of  a  unidirectional  lamina  can  be  expressed  as  (Dvorak,  1992) 


AE  +  X  DnP*> 

(7) 

s=hQ 

Br°+  X  FrA>  r=\,Q. 

(8) 

s=\,Q 


Here,  Frs  and  D  represent  stress  and  the  strain  influence  functions. 

The  phase  transfonnation  fields,  A,p ,  appearing  in  Eqs.  (2),  (6-8)  are  generally  caused  by 
local  phenomena  in  response  to  other  induced  fields.  Considering  a  piezo/pyro-electrically  active 
phase  subjected  to  electrical  field  <E  and  temperature  change  6 ,  the  transformation  stresses  and 
strains  caused  in  this  phase  can  be  written  as  , 

A--eT<E-Lad,  ju  =  dTE  +  ad,  (9) 
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where  a  is  the  coefficient  of  thermal  expansion,  and  e,d ,  denote  piezoelectric  constants.  For 
transverse  isotropy  with  x3  the  axis  of  rotational,  they  take  the  following  fonn  (Bahei-El-Din, 
2009;  Chen,  1994) 


"  0 

0 

d3l 

~  0 

0 

en 

aT 

0 

0 

d3\ 

0 

0 

en 

aT 

0 

0 

^15 

d33 

,  eJ  = 

0 

0 

en 

,  a  =  < 

aL 

0 

0 

0 

0 

0 

^15 

0 

0 

ei5 

0 

0 

0 

0 

0 

0  _ 

0 

0 

0 

0 

(10) 


where  aL  and  aT  denote  longitudinal  and  transverse  coefficient  of  thermal  expansion, 
respectively.  The  piezoelectric  constants  are  related  by  (Bahei-El-Din,  2009) 

e31  =  2  kd3l  +  ld33,  e33  =  2  td3l  +  nd33,  e15  =  pdl5.  (11) 


The  formulation  is  completed  by  listing  the  electric  displacement  <D  and  electrical  field  <E 
as  the  superposition  of  the  direct  effect  and  the  coupling  effects, 


(D  =  k£E  +  ee-q9  =  kJE  +  d<j-qO, 

(12) 

<E  =  -  gcr  +  HO, 

(13) 

such  that, 

g  =  K~ld,  n  =  K*lq  (14) 


Here,  g  is  the  voltage  constant,  kct  and  k,  are  (3x3)  permittivity  matrices  measured  at  constant 
stress  or  constant  strain,  respectively,  and  q  is  pyroelectric  constant.  The  piezoelectric  constants 
are  inter-related  by  (Bahei-El-Din,  2009) 

Bn  =  dn /k33ct,  #33=  d33!  k33ct,  gl5  =  dl5 /kU(T.  (15) 


For  a  transversely  isotropic  phase,  the  permittivity  and  the  pyroelectric  constant  take  the  form 
(Chen,  1994;  Kumar,  2009) 


K11 

0 

0  ' 

o' 

K  = 

0 

Ku 

0 

»  cl^< 

0 

0 

0 

K33_ 

43. 

(16) 
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2.2.  MACRO-SCALE  MECHANICS  OF  A  FIBROUS  LAMINATE 


The  laminate  considered  consists  of  n  unidirectionally  reinforced  fibrous  composite 
laminas  with  a  total  thickness  h  and  lamina  thickness  1 1  ,  i  =1  ,..n  .  Hence,  ^ t!  =  h  and  the 

i  =1  ,n 

lamina  volume  fraction  is  cj  =  t.  /  h  .  The  global  system  of  coordinates  of  the  laminate  is  denoted 
by  x  j ,  j  =  1,2,3,  where  xtx2  is  mid-plane  of  the  laminate  and  x3  is  the  out  of  plane  direction. 
The  laminas  local  coordinate  system  is  denoted  by  xk  ,  k  =  1,2,3,  where  the  fibers  are  aligned  in 
the  x"3  direction  as  shown  in  Fig.  1.  The  fiber  orientation  is  specified  by  the  angle,  cp  ,  between 
the  x3  and  x1  axes.  The  laminate  is  subjected  to  membrane  forces  A  =  [Au,  A22,  A12]  and 
bending  moments  A  =  [An,  A22,  A12] .  In  addition,  the  laminas  may  be  subjected  to  a  change  of 

temperature  0 .  On  the  phase  scale,  fibers  may  possess  a  three-way,  thermo-electro-mechanical 
coupling  in  response  to  the  above  mentioned  mechanical  loads  and  thennal  change  in  addition  to 
an  electric  field  intensity  <£  =  (<£1,‘£2,‘£3) . 

The  membrane  forces  A  =  [An, A22, A12]  and  bending  moments  A  =  [Au,  A22,  A12] 

applied  to  a  fibrous  laminate  (Fig.  1)  can  be  related  to  the  stresses  caused  in  the  individual 
laminas  as  follows  (Bahei-El-Din  et  al.,  2010), 


hi  2 

A  =  [  &(z  )dz 

-A/2 

h  12 

M=  J  z  <t(z  )dz 

-h!  2 


=  z 


A,  -M,./2  ^ 

|  (j(z  )dz  , 

\zi-ti/2  y 

fr'+f/2  .  1 

!  I  z  cr(z  )dz 


<=l.n  \zt-t,/2 


J 


(17) 


(18) 


Equations  (17,18)  are  written  with  reference  to  the  global  coordinate  system  x  y ,  /  =1,2,3 
shown  in  Fig.  1,  where  z  refers  to  the  x3  coordinate  of  the  midplane  of  a  lamina.  Symbols 
decorated  with  a  top  hat  denote  a  list  of  in-plane  components.  Hence,  the  lamina  stress  vector 
appearing  in  Eqs.  (17,  18)  is  written  as  &  =  [an,a22,al2]- 


In-plane  strains  i(z)  =  [en,e2 2,2en]  are  assumed  to  vary  linearly  across  the  laminate 
thickness  (Bahei-El-Din  et  al.,  2010).  Hence, 

s(z)  =  e0+ztc ,  (19) 


where  e0  =[£°l,£°2,2s°2]  and  k  =  [/y  ,/c2,/y2]  are  the  strain  and  curvature  at  the  mid-plane  of  the 
laminate,  respectively.  They  are  related  to  the  applied  loads  by  (Bahei-El-Din  et  al.,  2010): 

£o=/rN  +  g'fM+f',  (20) 


7 


K  =  C'N  +  V'M  +  g'. 

(21) 

The  first  two  terms  in  Eqs.  (20,  21)  represent  the  contribution  of  the  mechanical  loads.  The  third 
tenn  represents  the  effects  caused  by  the  lamina  eignstresses. 

Matrices  A'  ,B,C\V  J'  ,g'  are  a  function  of  the  elastic  properties  of  the  laminas,  their 
volume  fractions  and  laminate  layup  (Bahei-El-Din  et  ah,  2010); 

A'  -  (I  -B'B)A~l,  B’  =  -A  BD' , 

(22) 

C  =  -VBA~\  V  =  \V -BA  XB\\ 

(23) 

f'  =  -B'g-A'f ,  g'  =  -C'f  -V'g, 

(24) 

where  A,B ,  £9,  /  ,  g  are  related  to  lamina  properties  by; 

i  =1  ,n  i  =1  ,n 

(25) 

V=  Z^i^2+zfV-’ 

i=l,n  / 

(26) 

ii 

M 

> 

Qrq 

II 

M 

> 

(27) 

i  =1  ,n  i  =1  ,n 


The  (3x3)  lamina  stiffness  matrix  L  can  be  expressed  in  terms  of  Hill’s  overall  elastic 
moduli  of  a  unidirectional  composite  as  follows  (Bahei-El-Din,  1992), 


L 


4km  2m  t  0 

2  mi  ELk+mn  0 

0  0  p(k  +m ) 


(28) 


In  analogy  with  Eq.  (8),  lamina  stresses  in  the  overall  coordinates  x  .  ,  are  written  as  the 

superposition  of  the  stresses  caused  by  the  laminate  loads  and  those  caused  by  the  lamina 
transformation  stresses; 

<T,=fiw+e,*+IiyV  (29) 

j  =1,« 


Here,  Pt ,  Qf  denote  stress  distribution  factors,  and  U y  denotes  stress  transformation  influence 

functions.  They  are  a  function  of  elastic  moduli  and  volume  fractions  of  the  laminas,  and  the 
laminate  layup  (Bahei-El-Din  et  ah,  2010); 
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Pi=Li(A'  +  ziC),  Qi=Li(0'+z  ,£>'), 


(30) 


Uu 


=  5uI-tjpi 


)Qr 


(31) 


Equilibrium  of  lamina  forces  leads  to  the  following  relations  (Bahei-El-Din  et  ah,  2010): 


Xmw,  'Z(t,zi)pi=o, 

i=\,n  i=\,n 


(32) 


1.1,0,  =«,  l(t,z,)e,=r- 

i  =l,n  i=\,n 


(33) 


For  completeness,  coordinate  transformation  for  the  lamina  stresses  and  strains  are  listed  as 
follows  (Bahei-El-Din  et  ah,  2010): 


=  R,<ri, 

At 

II 

,  Ai  = 

R,At,  Hi 

L, 

=  NfLiNi,  M. 

=  R-M1 

R,, 

cos2  <pt 

sin2  (pt 

-  \  sin  2<pi 

R- 

=  N7'  = 

sin2  <pj 

cos2  <pt 

2  sin  2cpi 

sin  2^;. 

-sin2^(. 

cos  2  <pt 

(34) 


(35) 


where  vectors  decorated  with  a  top  bar  are  referred  to  the  lamina  material  axes,  otherwise 
reference  is  made  to  the  laminate  overall  coordinates,  and  (pi  denotes  the  fiber  orientation,  Fig.  1. 

2.3.  GLOBAL-SCALE  ANALYSIS  OF  FIBROUS  COMPOSITE  STRUCTURES 

In  utilizing  the  transformation  field  analysis  for  composites  in  explicit  formulation  of  the 
finite  element  method,  it  is  essential  to  note  the  analogy  between  treatment  of  the  latter  for  non¬ 
mechanical  deformations,  e.g.  thermal  strains,  and  the  TFA  formulation.  In  both  cases,  the 
structure  is  treated  as  elastic,  with  superimposed  defonnations  stemming  from  the  non¬ 
mechanical  fields.  In  the  present  problem  for  piezoelectric  composites,  the  non-mechanical  strain 
and  stresses  are  due  to  electromechanical  coupling  of  the  phases.  Hence,  in  implementing  the 
TFA  analysis  into  any  finite  element  framework,  the  eigen  stresses  are  computed  with  the  above 
formulation  and  its  underlying  micromechanical  model  and  used  to  update  the  nodal  loads. 

For  an  eigen  stress  Am  found  in  a  laminated  composite  element  of  a  finite  element  mesh, 
which  represents  some  structure,  the  nodal  forces  qm  ,  which  equilibrate  the  eigen  stresses  are 
given  by  the  principle  of  virtual  work, 
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(36) 


q  =  f  r  A  dV 

l  m  Jy  mm 


Where  r m  represents  spatial  derivatives  of  the  displacement  field  in  element  in  ,  and  V m  is 
volume  of  the  element. 

The  entire  TFA  model  is  coded  into  the  user  defined  subroutine  UEXPAN  of  the  ABAQUS 
finite  element  program.  This  serves  as  the  interface  between  the  global  finite  element  mesh  and 
the  response  of  the  individual  elements.  The  latter  are  symmetric  laminates  for  which  the  eigen 
strains  and  stresses  are  calculated  with  the  TFA  method  described  in  Sections  2.1  and  2.2  and 
supplied  to  the  finite  element  scheme. 


3.  RESULTS 

In  preparation  for  a  finite  element  analysis  of  a  laminated  structure,  the  TFA  scheme  is  used 
first  to  prepare  the  data  required  in  the  finite  element  solution,  namely  the  overall  elastic 
properties  of  the  laminate,  and  the  data  required  for  the  transformation  field  analysis,  namely  the 
stress  and  strain  concentration  factors  and  influence  functions.  In  this  pre-processing  stage,  a 
micromechanical  model  is  to  be  utilized  for  the  micromechanical  formulation  of  Section  2.1.  In 
the  present  work,  a  refined,  periodic  hexagonal  array  (PHA)  model  (Dvorak  and  Teply,  1985)  is 
applied.  In  this  case,  a  representative  volume  element  (RVE)  is  divided  into  subvolumes  F  , 
r  =  \,Q,  and  analyzed  under  certain  loads,  which  are  derived  from  eqs.  (7,  8)  to  compute  the 
columns  of  the  concentration  factors  Ar ,  B  ,  and  influence  functions  D  ,  Frs ,  r,s  =  l,Q  . 

Two  applications  are  presented;  one  is  a  benchmark,  classical  laminated  plate  problem,  and 
one  which  involves  a  laminated  plate  with  complex  geometry.  In  both  problems  the 
unidirectional  composite  used  consists  of  piezoelectric  fibers  (PZT-5A)  embedded  in  an  epoxy 
(DY-063)  matrix.  The  mechanical  properties  of  both  constituents  are  given  in  Table  1,  and 
electrical  properties  of  the  fibers  are  given  in  Table  2.  The  fiber  volume  fraction  is  kept  at  0.55 
for  all  plies.  Using  the  laminated  plate  analysis  (Section  2.2)  with  the  unidirectional  plies 
represented  by  the  PHA  model  and  an  RVE  with  60  fiber  subvolumes  and  240  matrix 
subvolumes,  the  overall  elastic  properties  are  determined  for  two  symmetric  laminates,  a  cross 
ply,  (0/90)s,  and  a  quasi-isotropic,  (0/±45/90)s  laminate.  The  results  are  given  in  Table  3. 

3.1.  BENCHMARK  PROBLEMS 

To  verily  that  the  tri-lengthscale  analysis  using  TFA/ ABAQUS  yields  the  same  results  as  the 
stand  alone  TFA,  a  solid  plate  constructed  from  a  (0/90)s  laminate  is  considered.  The  electrically 
active  fibers  in  the  0°  lamina  are  subjected  to  an  electric  field  cEi  in  the  axial  direction  at  a 

magnitude  of  l.OxlO6  V/m.  This  induces  eigen  stresses  Xu  =  An  =  -e31<E3,  A33  =  -e33(E3  in  the  0°, 
PZT  fiber,  and  the  corresponding  overall  strain  at  mid-plane  of  the  laminate  is  computed  with  the 
TFA  scheme  as  =  2.70x10 4 ,  £•„  = -5.1 1x10  ? ,  and  s,2  =  0  .  Using  the  ABAQUS  finite 

element  code,  the  same  plate  problem  was  solved,  where  the  non-mechanical,  eigen  strains 
caused  by  the  electric  field  in  the  0°,  PZT  fiber  are  computed  through  the  user  subroutine 
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UEXPAN  with  the  TFA  scheme.  The  plate,  which  measured  0. 5x0.5  m,  was  subdivided  into  25 
quadrilateral  elements  of  equal  dimensions  and  was  unconstrained  at  the  boundaries.  The  finite 
element  results  obtained  in  this  case  were  identical  with  those  provided  by  the  laminate  analysis. 

3.2.  EXAMPLES 

Consider  a  square,  quasi-isotropic,  (0/±45/90)s  laminated  plate  with  a  concentric  circular 
hole.  The  side  length  of  the  plate  is  0.5  m,  and  diameter  of  the  hole  is  0.1  m.  The  tri-lengthscale 
TFA/ABAQUS  model  was  used  to  solve  the  problem,  where  the  underlying  micromechanical 
model  is  the  periodic  hexagonal  array.  As  indicated  earlier,  the  overall  elastic  properties  of  the 
laminate  together  with  the  stress  concentration  factors  and  influence  functions  of  the  fiber  and 
matrix  subvolumes  are  pre-determined  and  stored  for  processing  by  the  user-defined  subroutine 
UEXPAN.  The  overall  elastic  properties  of  the  quasi-isotropic  laminate  under  consideration  are 
given  in  Table  3. 

In  the  finite  element  analysis  one  quarter  of  the  plate  was  considered  with  boundary 
conditions  which  simulate  symmetry  of  the  overall  deformations.  In  one  problem,  the  plate  is 
otherwise  unconstrained  and  subjected  to  an  overall  tensile  stress  of  50  MPa  in  direction  of  the  0° 
fiber.  The  stress  contours  obtained  from  ABAQUS  for  an,  the  stress  in  direction  of  the  load  are 
shown  in  Fig.  2.  In  a  post-processor  exercise,  the  TFA  approach  was  used  to  determine  the  local 
fields  in  the  fiber  and  matrix  within  all  laminas  using  the  overall  stresses  found  from  the  finite 
element  analysis.  In  particular,  the  electric  displacement  generated  in  the  PZT  fiber  within  the  0° 
lamina  are  of  interest  as  it  reflects  the  stress  concentration  caused  by  discontinuity  of  the  plate 
geometry.  This  is  indicated  in  Fig.  3,  where  values  of  the  electric  displacement  at  selected 
locations  are  shown  as  insets  on  the  map  for  the  overall  stress  in  direction  of  the  applied  load.  It 
is  seen  that  the  magnitude  of  the  computed  electric  displacement  reflects  the  stress  gradient.  It  is 
worth  mentioning  that  the  electric  displacement  in  the  same  laminated  plate  without  the  circular 
discontinuity  at  overall  stress  of  50  MPa  is  6.3x10'  C/nr.  The  same  magnitude  is  found  in  the 
plate  under  consideration  at  the  upper-right  comer  where  the  stress  field  is  closest  to  that  found 
in  the  solid  plate.  On  the  other  hand,  the  electric  displacement  found  at  perimeter  of  the  circular 
cutout  is  2.2  C/m”  at  the  stress  free,  three-o’clock  location.  This  is  caused  by  the  overall  stress  in 
the  X2  direction,  which  is  transverse  to  the  load.  The  electric  displacement  at  the  twelve-o’clock 
location  where  the  stress  concentration  is  the  largest  is  18.8  C/m“. 

In  another  loading  condition,  the  PZT  fiber  in  the  0°  plies  is  subjected  to  an  electric  field  E3 

of  l.OxlO6  V/m  in  the  longitudinal  direction.  In  this  problem,  the  boundary  conditions  are 
modified  to  impose  a  displacement  constraint  in  the  direction  perpendicular  to  the  outer  edges  of 
the  plate  in  addition  to  the  symmetry  boundary  conditions.  Following  a  finite  element  analysis 
with  the  ABAQUS  program  with  the  TFA  interface,  post  processing  of  the  results  provides  the 
ply  and  phase  stresses.  To  illustrate  this  capability  of  the  presented  tri-multiscale  analysis,  we  list 


11 


in  Table  4  the  ply  stresses,  <j1i  ,  a22 ,  a  12  ,  and  the  fiber  axial  stress,  erf,  at  locations  where  each 
of  the  overall  stresses  Su,  S22,  S12  is  a  maximum. 

Additional  post-processing  of  the  finite  element  results  using  the  micromechanical/laminate 
analysis  of  Sections  2.1  and  2.2  yields  the  electric  displacement  in  the  active  fibers  of  the  0°  plies 
caused  by  the  applied  electrical  field,  <E} ,  introduced  in  the  same  fiber.  The  electric  displacement 

thus  calculated  is  due  to  the  overall  stresses  generated  as  a  result  of  the  plate  geometry  and 
constraints.  When  this  is  superimposed  on  the  electric  displacement  due  to  (1)  the  direct  electric 
effect,  k33<E3  ,  and  (2)  the  local  stress  generated  as  a  result  of  the  mutual  constraints  of  the  phases 

and  plies,  and  the  result  is  divided  by  the  applied  electric  field,  the  apparent  pennittivity  is 
obtained.  Figure  4  shows  the  apparent  permittivity  at  selected  locations,  superimposed  on 
contours  of  the  overall  stress  found  in  direction  of  the  Xi  axis.  We  note  that  the  pennittivity  is 
the  smallest  at  the  twelve-o’clock  location  at  perimeter  of  the  circular  hole  where  the  stress  due 
to  the  electric  field  in  the  0°  fiber  is  the  largest.  In  contrast,  the  permittivity  at  the  three-o’clock 
location  is  9.6x1 0‘9  C/V.m.  As  expected,  this  is  close  to  that  found  in  a  laminated  plate  of  the 
same  geometry  and  material  (Bahei-El-Din  and  Micheal,  2012). 


4.  CONCLUSIONS 

A  tri-lengthscale  model  for  analysis  of  composite  structures  with  piezoelectric  constituents 
is  presented.  It  combines  a  micromechanical  model  for  fibrous  composites,  a  laminated  plate 
analysis  under  plane  stress,  and  a  finite  element  analysis.  The  micromechanical  model  utilized  is 
a  refined  one  in  which  the  stress  and  strain  fields  in  the  fiber  and  matrix  are  represented  by 
piecewise  unifonn  idealizations.  The  finite  element  analysis  was  perfonned  with  the  ABAQUS 
program.  Non-mechanical  stresses  and  strains  generated  in  composite  phases  are  treated  as 
transfonnation  (eigen)  fields  and  their  effect  on  the  composite  media  is  determined  with  a 
transformation  field  analysis. 

The  applications  presented  cover  a  benchmark  problem,  which  verifies  that  the  results  of  the 
tri-lengthscale  TF AJ ABAQUS  model  are  reliable,  and  a  plate  with  a  hole,  which  illustrates  the 
capabilities  of  the  solution  scheme  in  analyzing  complex  geometries  with  detailed  underlying 
models.  Extension  of  this  solution  strategy  to  other  eigen  fields,  stemming  for  example  from 
pyroelectric  effects  and  damage,  is  quite  feasible  and  currently  in  progress  by  the  authors. 
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Table  1.  Mechanical  properties  of  the  fiber  and  matrix  (Bahei-El-Din,  2009;  Berlincourt  et 
ah,  2009). 


Material 

E  (GPa) 

V 

PZT-5A 

60 

0.34 

DY063  Epoxy 

3.35 

0.35 

Table  2. 

Piezo/pyroelectric  properties  of  PZT-5A  fibers  (Bahei-El-Din,  2009;  Berlincourt 
et  ah,  2009). 

dj,  i 

J33 

d$  i 

3 31 

$33 

15 

(10‘12  m/V) 

(10‘3  V.m/N) 

-171 

347 

584 

11.4 

24.8 

38.2 

Table  3.  Computed  overall  mechanical  properties  of  laminates. 

Layup  En  (GPa) 

E22  (GPa) 

G12  (GPa) 

V12 

(0/90)s  22.37 

22.37 

3.69 

0.15 

(0/±45/90)s  17.76 

17.76 

6.71 

0.325 
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Table  4.  Lamina  and  fiber  stresses  at  locations  of  maximum  overall  stress  components 
(MPa). 


Overall 

Stress 

Lamina 

CT// 

^22 

°12 

0° 

-14.60 

-.21 

0.48 

-25.4 

45° 

-5.12 

-1.10 

-0.64 

-6.02 

s„ 

-45° 

-6.79 

-2.77 

-3.18 

-13.5 

vo 

o 

o 

-3.86 

-2.91 

0.48 

5.87 

0° 

-9.58 

-.57 

1.08 

-16.5 

45° 

-2.38 

-.11 

1.35 

0.72 

s,2 

-45° 

-6.12 

-3.85 

4.37 

-16.1 

VO 

o 

o 

-2.64 

.33 

1.08 

1.12 

0° 

-1.06 

1.28 

0.02 

-2.10 

45° 

0.65 

2.01 

0.69 

3.36 

s22 

-45° 

0.57 

1.92 

-0.56 

3.0 

vo 

o 

o 

0.05 

4.87 

0.02 

8.46 
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Fig.  1.  Geometry  and  load  of  a  fibrous  laminate. 
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PZT-5A/DY063,  (0/+45/90)s  ,  C,=  0.55 


Fig.  2.  Contours  of  an  stress  (MPa)  caused  by  an  overall  stress  of 
50  MPa  in  the  Xi  direction. 
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PZT-5A/DY063,  (0/±45/90)s  ,  c,=  0.55 


Fig.  3.  Electric  displacement  <D3  computed  in  0°  PZT  fiber  (10"2  C/m2) 

superimposed  on  an  stress  contours  (Fig.  2)  at  selected  locations, 
due  to  an  overall  stress  of  50  MPa  in  the  Xi  direction. 
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PZT-5A/DY063,  (0/±45/90)s,  Cf  =  0.55 


s,  Sll 


(Av«.  Crit. :  75%) 


I 


+7.22e-02 
-5. 67e-01 
-1.  21o  +  00 
-1.  85e  +  00 
-2. 48e+00 
-3. 12e+00 
-3.  76e  +  00 
-4. 4  0e  +  00 
-5. 04e+00 
-5. 68e+00 
-6.  32e+00 
-6.  96e  +  00 
-7. 60e+00 


Fig.  4.  Apparent  permittivity  computed  in  0°  PZT  fiber  (10"2  C/V.m) 
superimposed  on  an  stress  contours  (MPa)  at  selected  locations, 
due  to  electric  field  cEi  of  l.OxlO6  V/m  induced  in  the  0°  fiber. 
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